# Seurat
library(Seurat)
# Data
library(dplyr)
# Plotting
library(ggplot2)
library(RColorBrewer)
library(patchwork)
library(gridExtra)
library(grid)
library(viridis)
library(ComplexHeatmap)
Attaching SeuratObject
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Attaching package: ‘gridExtra’
The following object is masked from ‘package:dplyr’:
combine
Loading required package: viridisLite
========================================
ComplexHeatmap version 2.8.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
genomic data. Bioinformatics 2016.
The new InteractiveComplexHeatmap package can directly export static
complex heatmaps into an interactive Shiny app with zero effort. Have a try!
This message can be suppressed by:
suppressPackageStartupMessages(library(ComplexHeatmap))
========================================
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
# Source files
source("plotting_global.R")
# Seurat files
so_file <- "data/object/seurat.rds"
so_out_file <- "data/object/seurat_clean.rds"
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()
so <- readRDS(so_file)
# Create clean Seurat object
so <- CreateSeuratObject(counts = GetAssayData(so, assay = "RNA", slot = "counts"), meta.data = dplyr::select(so@meta.data, -rna_snn_res.0.8))
so_l <- SplitObject(so, split.by = "treatment")
so_l <- lapply(so_l, NormalizeData)
so_l <- lapply(so_l, FindVariableFeatures, nfeatures = 2000)
cc_kowalczyk <- read.csv("cc_kowalczyk.csv")
cc_kowalczyk <- cc_kowalczyk[cc_kowalczyk$sig_population >= 5, ]
variable_features_reduce <- function(so, genes) {
variable_features <- VariableFeatures(so)
VariableFeatures(so) <- variable_features[!variable_features %in% genes]
return(so)
}
so_l <- lapply(so_l, variable_features_reduce, genes = cc_kowalczyk$gene)
so_l <- lapply(so_l, ScaleData, vars.to.regress = c("nCount_RNA"))
Regressing out nCount_RNA Centering and scaling data matrix Regressing out nCount_RNA Centering and scaling data matrix
options(warn = -1)
so_l <- lapply(so_l, RunPCA, npcs = 100, verbose = FALSE)
so_l <- lapply(so_l, FindNeighbors, dims = 1:30, verbose = FALSE)
so_l <- lapply(so_l, FindClusters, verbose = FALSE)
so_l <- lapply(so_l, RunUMAP, dims = 1:90, verbose = FALSE)
so_l <- lapply(so_l, RunTSNE, dims = 1:30, verbose = FALSE)
reduction <- "umap"
options(repr.plot.width = 25, repr.plot.height = 20)
plot <- function(so) {
dplot_1 <- DimPlot(so, reduction = reduction, group.by = "RNA_snn_res.0.8", label = TRUE) &
ggtitle("Cluster") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "none") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_2 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) &
ggtitle("Tissue") & xlab("UMAP 1") & ylab("UMAP 2") &
scale_color_manual(values = color$tissue, na.value = "dark gray") &
theme(aspect.ratio = 1, legend.position = "bottom") &
guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
dplot_3 <- DimPlot(so, reduction = reduction, group.by = "sample_rep", label = FALSE) &
ggtitle("Replicate") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
dplot_4 <- DimPlot(so, reduction = reduction, group.by = "main_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_5 <- DimPlot(so, reduction = reduction, group.by = "fine_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
fplot_1 <- FeaturePlot(so, reduction = reduction, features = "Cd34") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so, reduction = reduction, features = "Kit") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so, reduction = reduction, features = "Gata2") & theme(aspect.ratio = 1)
fplot_4 <- FeaturePlot(so, reduction = reduction, features = "Gata1") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so, reduction = reduction, features = "Tfrc") & theme(aspect.ratio = 1) & ggtitle("Tfrc (CD71)")
fplot_6 <- FeaturePlot(so, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
fplot_7 <- FeaturePlot(so, reduction = reduction, features = "msTcr_RNA1") & theme(aspect.ratio = 1)
fplot_8 <- FeaturePlot(so, reduction = reduction, features = "msIgkc_RNA1") & theme(aspect.ratio = 1)
fplot_9 <- FeaturePlot(so, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
fplot_10 <- FeaturePlot(so, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
fplot_11 <- FeaturePlot(so, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)
dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + fplot_7 + fplot_8 + fplot_9 + fplot_10 + fplot_11 + plot_layout(ncol = 4)
}
lapply(so_l, plot)
$CpG $NaCl
so_cpg <- so_l[[1]]
so_nacl <- so_l[[2]]
CpG: Cluster 17
NaCl: Cluster 17
so_cpg <- subset(so_cpg, subset = RNA_snn_res.0.8 != 17)
so_nacl <- subset(so_nacl, subset = RNA_snn_res.0.8 != 17)
CpG: None
NaCl: Cluster 18
so_nacl <- subset(so_nacl, subset = RNA_snn_res.0.8 != 18)
options(repr.plot.width = 10, repr.plot.height = 5)
tissue_cpg_cluster <- ggplot(so_cpg@meta.data, aes(x = seurat_clusters, fill = tissue)) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$tissue) +
ggtitle("CpG FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom")
tissue_nacl_cluster <- ggplot(so_nacl@meta.data, aes(x = seurat_clusters, fill = tissue)) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$tissue) +
ggtitle("NaCl FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom")
tissue_cpg_cluster + tissue_nacl_cluster + plot_layout(ncol = 2)
clov: cluster overlap
clpu: cluster pure
clmi: cluster mixed
cpg_clov <- dplyr::group_by(so_cpg@meta.data, RNA_snn_res.0.8) %>% count(tissue)
# Get all pure cluster with no FACS mixture
cpg_clpu <- cpg_clov[!cpg_clov$RNA_snn_res.0.8 %in% cpg_clov$RNA_snn_res.0.8[duplicated(cpg_clov$RNA_snn_res.0.8)], ]
# Get all mixed cluster
cpg_clmi <- cpg_clov[!cpg_clov$RNA_snn_res.0.8 %in% cpg_clpu$RNA_snn_res.0.8, ]
# Remove cells with low frequncy mixed cluster from Seurat object
cpg_clmi_lf <- cpg_clmi[cpg_clmi$n <= 2, ] # Get all low frequency mixed cluster
so_cpg$clmi_lf <- ifelse(paste0(so_cpg$`RNA_snn_res.0.8`, so_cpg$tissue) %in% paste0(cpg_clmi_lf$RNA_snn_res.0.8, cpg_clmi_lf$tissue), TRUE, FALSE)
so_cpg <- subset(so_cpg, subset = clmi_lf == FALSE)
# Get all mixed cluster with high frequency for myeloid dominant and progenitor dominant
cpg_clmi_hf <- cpg_clmi[!cpg_clmi$RNA_snn_res.0.8 %in% cpg_clmi_lf$RNA_snn_res.0.8, ]
cpg_clmi_hf <- dplyr::group_by(cpg_clmi_hf, RNA_snn_res.0.8) %>% dplyr::filter(n == min(n))
cpg_clmi_hf_m <- cpg_clmi_hf[cpg_clmi_hf$tissue == "Myeloid", ]
cpg_clmi_hf_p <- cpg_clmi_hf[cpg_clmi_hf$tissue == "Progenitor", ]
nacl_clov <- dplyr::group_by(so_nacl@meta.data, RNA_snn_res.0.8) %>% count(tissue)
# Get all pure cluster with no FACS mixture
nacl_clpu <- nacl_clov[!nacl_clov$RNA_snn_res.0.8 %in% nacl_clov$RNA_snn_res.0.8[duplicated(nacl_clov$RNA_snn_res.0.8)], ]
# Get all mixed cluster
nacl_clmi <- nacl_clov[!nacl_clov$RNA_snn_res.0.8 %in% nacl_clpu$RNA_snn_res.0.8, ]
# Remove cells with low frequncy mixed cluster from Seurat object
nacl_clmi_lf <- nacl_clmi[nacl_clmi$n <= 2, ] # Get all low frequency mixed cluster
so_nacl$clmi_lf <- ifelse(paste0(so_nacl$`RNA_snn_res.0.8`, so_nacl$tissue) %in% paste0(nacl_clmi_lf$RNA_snn_res.0.8, nacl_clmi_lf$tissue), TRUE, FALSE)
so_nacl <- subset(so_nacl, subset = clmi_lf == FALSE)
# Get all mixed cluster with high frequency for myeloid dominant and progenitor dominant
nacl_clmi_hf <- nacl_clmi[!nacl_clmi$RNA_snn_res.0.8 %in% nacl_clmi_lf$RNA_snn_res.0.8, ]
nacl_clmi_hf <- dplyr::group_by(nacl_clmi_hf, RNA_snn_res.0.8) %>% dplyr::filter(n == min(n))
nacl_clmi_hf_m <- nacl_clmi_hf[nacl_clmi_hf$tissue == "Myeloid", ]
nacl_clmi_hf_p <- nacl_clmi_hf[nacl_clmi_hf$tissue == "Progenitor", ]
find_tissue_markers <- function(so, cluster, low_freq, high_freq, min_cells_group = 3) {
so <- subset(so, subset = RNA_snn_res.0.8 %in% cluster)
# Set idents to cluster tissue combination
so$cluster_tissue <- paste0(so$`RNA_snn_res.0.8`, "_", so$tissue)
Idents(so) <- "cluster_tissue"
# List to store results from looping through cluster
deg_l <- list()
for(cluster_idx in unique(so$`RNA_snn_res.0.8`)) {
# Check if enough cells are present in both idents
so_cluster_idx <- subset(so, RNA_snn_res.0.8 == cluster_idx)
sample_size_check <- all(table(so_cluster_idx$tissue) > min_cells_group) & length(unique(so_cluster_idx$tissue)) == 2
# DEG
if(sample_size_check) {
deg <- FindMarkers(so, ident.1 = paste0(cluster_idx, "_", low_freq), ident.2 = paste0(cluster_idx, "_", high_freq), min.cells.group = min_cells_group, only.pos = TRUE, verbose = TRUE)
# Filter deg list
if(nrow(deg) > 0) {
deg <- deg[grepl("^(mt-|Rp(s|l)|Gm\\d)", rownames(deg)) == FALSE, ]
deg <- deg[!rownames(deg) %in% cc_kowalczyk$gene, ]
}
deg_l[[cluster_idx]] <- deg
} else {
deg_l[[cluster_idx]] <- NULL
}
}
return(deg_l)
}
deg_cpg_clmi_hf_m <- find_tissue_markers(so_cpg, cpg_clmi_hf_m$RNA_snn_res.0.8, low_freq = "Myeloid", high_freq = "Progenitor")
deg_cpg_clmi_hf_p <- find_tissue_markers(so_cpg, cpg_clmi_hf_p$RNA_snn_res.0.8, low_freq = "Progenitor", high_freq = "Myeloid")
deg_nacl_clmi_hf_m <- find_tissue_markers(so_nacl, nacl_clmi_hf_m$RNA_snn_res.0.8, low_freq = "Myeloid", high_freq = "Progenitor")
# deg_nacl_clmi_hf_p <- find_tissue_markers(so_nacl, nacl_clmi_hf_p$RNA_snn_res.0.8, low_freq = "Progenitor", high_freq = "Myeloid")
top_hm <- function(cluster_idx, so, deg, top = 25) {
# Filter deg list
deg <- deg[[cluster_idx]]
# Get up and downreagulated genes
deg_up <- deg[sign(deg$avg_log2FC) == 1, ]
deg_down <- deg[sign(deg$avg_log2FC) == -1, ]
# Select top hits by adjusted p value
deg_up <- arrange(deg_up, p_val_adj)[1:top, ]
deg_down <- arrange(deg_down, p_val_adj)[1:top, ]
# Combine top hits
deg <- rbind(deg_up, deg_down)
# Extract normalized count matrix for top hits
so <- subset(so, subset = RNA_snn_res.0.8 == cluster_idx)
cnt_norm <- GetAssayData(so, assay = "RNA", slot = "data")
cnt_norm <- cnt_norm[rownames(cnt_norm) %in% rownames(deg), ]
# Order count matrix by tissue
meta <- so@meta.data
meta <- arrange(meta, tissue)
cnt_norm <- cnt_norm[, rownames(meta)]
# Compute mean per group
colnames(cnt_norm) <- meta$tissue
cnt_norm <- data.frame(
Myeloid = rowMeans(cnt_norm[, grepl("Myeloid", colnames(cnt_norm), fixed = TRUE)]),
Progenitor = rowMeans(cnt_norm[, grepl("Progenitor", colnames(cnt_norm), fixed = TRUE)]))
# ComplexHeatmap color ramp
color_range <- max(abs(cnt_norm))
color_break <- seq(0, color_range, 0.01)
color_ramp <- viridis(length(color_break), option = "magma")
# Top annotaion
annoation_col <- data.frame("FACS" = colnames(cnt_norm))
rownames(annoation_col) <- colnames(cnt_norm)
annotation_colors <- list("FACS" = color$tissue)
# Heat map
hm <- grid.grabExpr(draw(ComplexHeatmap::pheatmap(
mat = as.matrix(cnt_norm),
main = paste("Cluster", cluster_idx),
fontsize_row = 10,
scale = "none",
cluster_rows = TRUE,
cluster_cols = FALSE,
cellwidth = 12,
cellheight = 12,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "complete",
show_row_dend = FALSE,
show_rownames = TRUE,
show_colnames = FALSE,
annotation_col = annoation_col,
annotation_colors = annotation_colors,
color = color_ramp,
breaks = color_break,
border_color = NA,
use_raster = TRUE)))
return(hm)
}
# CpG
hm_cpg_clmi_hf_m <- lapply(names(deg_cpg_clmi_hf_m), top_hm, so = so_cpg, deg = deg_cpg_clmi_hf_m)
hm_cpg_clmi_hf_p <- lapply(names(deg_cpg_clmi_hf_p), top_hm, so = so_cpg, deg = deg_cpg_clmi_hf_p)
# NaCl
hm_nacl_clmi_hf_m <- lapply(names(deg_nacl_clmi_hf_m), top_hm, so = so_nacl, deg = deg_nacl_clmi_hf_m)
options(repr.plot.width = 20, repr.plot.height = 10)
gridExtra::arrangeGrob(grobs = hm_cpg_clmi_hf_m, nrow = 2) %>% grid::grid.draw()
so_cpg$clmi_hf <- ifelse(paste0(so_cpg$`RNA_snn_res.0.8`, so_cpg$tissue) %in% paste0(c(2, 8, 1, 4, 0, 6, 3, 12, 7), "Myeloid"), TRUE, FALSE)
so_cpg <- subset(so_cpg, subset = clmi_hf == FALSE)
Keep all. It looks like a mixture of myeloid and MEP progenitors
options(repr.plot.width = 20, repr.plot.height = 5)
gridExtra::arrangeGrob(grobs = hm_cpg_clmi_hf_p, nrow = 1) %>% grid::grid.draw()
Only keep cluster 14 because it looks like a true myeloid progenitor cluster
options(repr.plot.width = 20, repr.plot.height = 10)
gridExtra::arrangeGrob(grobs = hm_nacl_clmi_hf_m, nrow = 2) %>% grid::grid.draw()
so_nacl$clmi_hf <- ifelse(paste0(so_nacl$`RNA_snn_res.0.8`, so_nacl$tissue) %in% paste0(c(5, 7, 12, 6, 10, 0), "Myeloid"), TRUE, FALSE)
so_nacl <- subset(so_nacl, subset = clmi_hf == FALSE)
so_l <- list(so_cpg, so_nacl)
so_l <- lapply(so_l, FindVariableFeatures, nfeatures = 2000)
so_l <- lapply(so_l, variable_features_reduce, genes = cc_kowalczyk$gene)
so_l <- lapply(so_l, ScaleData, vars.to.regress = c("nCount_RNA"))
Regressing out nCount_RNA Centering and scaling data matrix Regressing out nCount_RNA Centering and scaling data matrix
options(warn = -1)
so_l <- lapply(so_l, RunPCA, npcs = 100, verbose = FALSE)
so_l <- lapply(so_l, FindNeighbors, dims = 1:30, verbose = FALSE)
so_l <- lapply(so_l, FindClusters, verbose = FALSE)
so_l <- lapply(so_l, RunUMAP, dims = 1:90, verbose = FALSE)
so_l <- lapply(so_l, RunTSNE, dims = 1:30, verbose = FALSE)
reduction <- "umap"
options(repr.plot.width = 25, repr.plot.height = 20)
plot <- function(so) {
dplot_1 <- DimPlot(so, reduction = reduction, group.by = "RNA_snn_res.0.8", label = TRUE) &
ggtitle("Cluster") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "none") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_2 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) &
ggtitle("Tissue") & xlab("UMAP 1") & ylab("UMAP 2") &
scale_color_manual(values = color$tissue, na.value = "dark gray") &
theme(aspect.ratio = 1, legend.position = "bottom") &
guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
dplot_3 <- DimPlot(so, reduction = reduction, group.by = "sample_rep", label = FALSE) &
ggtitle("Replicate") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
dplot_4 <- DimPlot(so, reduction = reduction, group.by = "main_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_5 <- DimPlot(so, reduction = reduction, group.by = "fine_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
fplot_1 <- FeaturePlot(so, reduction = reduction, features = "Cd34") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so, reduction = reduction, features = "Kit") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so, reduction = reduction, features = "Gata2") & theme(aspect.ratio = 1)
fplot_4 <- FeaturePlot(so, reduction = reduction, features = "Gata1") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so, reduction = reduction, features = "Tfrc") & theme(aspect.ratio = 1) & ggtitle("Tfrc (CD71)")
fplot_6 <- FeaturePlot(so, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
fplot_7 <- FeaturePlot(so, reduction = reduction, features = "msTcr_RNA1") & theme(aspect.ratio = 1)
fplot_8 <- FeaturePlot(so, reduction = reduction, features = "msIgkc_RNA1") & theme(aspect.ratio = 1)
fplot_9 <- FeaturePlot(so, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
fplot_10 <- FeaturePlot(so, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
fplot_11 <- FeaturePlot(so, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)
dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + fplot_7 + fplot_8 + fplot_9 + fplot_10 + fplot_11 + plot_layout(ncol = 4)
}
lapply(so_l, plot)
[[1]] [[2]]
options(repr.plot.width = 10, repr.plot.height = 5)
tissue_cpg_cluster <- ggplot(so_l[[1]]@meta.data, aes(x = seurat_clusters, fill = tissue)) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$tissue) +
ggtitle("CpG FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom")
tissue_nacl_cluster <- ggplot(so_l[[2]]@meta.data, aes(x = seurat_clusters, fill = tissue)) +
geom_bar(stat = "count", position = "fill") +
scale_fill_manual(values = color$tissue) +
ggtitle("NaCl FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") +
theme(legend.position = "bottom")
tissue_cpg_cluster + tissue_nacl_cluster + plot_layout(ncol = 2)
so <- merge(so_l[[1]], so_l[[2]])
saveRDS(so, so_out_file)
so <- NormalizeData(
object = so,
assay = "RNA",
normalization.method = "LogNormalize",
scale.factor = 10000
)
so <- FindVariableFeatures(
object = so,
assay = "RNA",
selection.method = "vst",
nfeatures = 3000
)
variable_features <- VariableFeatures(so)
VariableFeatures(so) <- variable_features[!variable_features %in% c(cc_kowalczyk$gene)]
so <- ScaleData(
object = so,
assay = "RNA",
vars.to.regress = c("nCount_RNA"),
do.scale = TRUE,
do.center = TRUE
)
Regressing out nCount_RNA Centering and scaling data matrix
source("bin/seurat_function.R")
so <- dim_clust(
# Seurat
so = so,
assay = "RNA",
# PCA features
blacklist_genes = NULL,
# Dim reduction
dims_pca = 100, # Default 20 - RunPCA dims
dims_umap = 90, # Default NULL - RunUMAP dims
dims_tsne = 30, # Default 1:5
min_dist = 0.3, # Default 0.3 - RunUMAP dmin.ist - controls how tightly the embedding
# Cluster
dims_cluster = 10, # Default 1:10 - FindNeighbors dims
cluster_res = 0.8 # Default 0.8 - FindClusters resoluton - above(below) 1.0 to obtain larger(smaller) number of communities)
)
PC_ 1 Positive: Car2, Prdx2, Blvrb, Glrx5, Rps2, Hbb-bs, Hba-a1, Cd24a, Hbb-bt, Cpox Ctse, Rhd, Car1, Mgst3, Hba-a2, Cldn13, Odc1, Ybx3, Ermap, Stard10 Hemgn, Gypa, C1qtnf12, Klf1, Rpl14, Hdgf, Hsp90aa1, Aqp1, Ubac1, Smim1 Negative: Tyrobp, Tmsb4x, Ptprc, Laptm5, Cd52, Fcer1g, Cyba, Ctss, Coro1a, Cd74 Lcp1, Lsp1, Cst3, Sh3bgrl3, Napsa, Arpc1b, Gm2a, Spi1, Lst1, Gngt2 Rgs2, Sat1, Psap, Fyb, Cybb, H2-Aa, Pou2f2, Pld4, Actb, H2-Ab1 PC_ 2 Positive: Fth1, Ftl1, Hba-a2, Hbb-bt, Hba-a1, Hbb-bs, Ctsb, Gypa, Hebp1, Slc4a1 Mgst3, Ctse, Hemgn, Blvrb, Cd24a, Rhd, Slc25a37, Cldn13, Hagh, Acp5 Cpox, Tfrc, Alas2, Ubac1, Prdx2, Isg20, C1qc, Ctsd, Vcam1, C1qa Negative: Myb, S100a10, Npm1, Tagln2, Gm15915, Ptprcap, Ldha, Cmtm7, Mif, Srm H2afy, Sox4, Phgdh, Nop10, Ybx1, Npm3, Actg1, Tmsb10, Ppp1r14b, Hsp90ab1 Mt1, Prdx6, Fbl, Hspd1, Hspa8, Srgn, Sptssa, Tomm5, Slc25a4, Psmb8 PC_ 3 Positive: Vcam1, C1qa, C1qc, C1qb, Fcna, Cd5l, Cadm1, Mertk, Igf1, Mrc1 Axl, Frmd4b, Itgad, Maf, Actn1, Id2, Slc40a1, Spic, Ccr3, Dnase1l3 Mafb, Hpgds, Asb2, Ptgs1, Cfp, Kcna2, Epb41l3, Slpi, Hmox1, Aif1 Negative: Ace, Cebpb, Itgal, Msrb1, Lyz2, Cyp4f18, Cx3cr1, Spn, Ceacam1, Cd300ld Tppp3, Smpdl3b, Nupr1, Hp, Klf4, Pglyrp1, Smpdl3a, Nr4a1, Eno3, Rbpms Fabp4, Ear2, Stk10, Ifitm2, Klf2, Tnfrsf1b, S1pr5, Il17ra, Rras, Dusp16 PC_ 4 Positive: Vcam1, C1qc, C1qa, Fcna, Mertk, Cd5l, C1qb, Igf1, Axl, Mafb Mrc1, Maf, Itgad, Slc40a1, Hmox1, Ccr3, Apoe, Slpi, Spic, Ptgs1 Epb41l3, Trf, Selenop, Kcna2, Kcnj10, Pla2g15, Clec4n, Sema6d, Timp2, Gfra2 Negative: Ciita, Ccnd1, H2-Eb1, H2-Oa, Dnase1l3, H2-Aa, H2-Ab1, Slamf7, H2-DMb2, Ppp1r14a Traf1, Lsp1, Tbc1d4, S100a11, Ffar2, Siglecg, Wdfy4, Jaml, Adam23, Cd74 Avpi1, Rogdi, Lyst, Flt3, Ccl5, S100a4, H2-Ob, Cd52, Gsn, Cd7 PC_ 5 Positive: Xcr1, Naaa, Rab7b, Irf8, Cd8a, Clec9a, Ppt1, P2ry14, Tlr3, Cadm1 A530099J19Rik, AC163354.1, Ifi205, Mpeg1, Pianp, Pdia5, Itgae, Gcsam, Olfm1, Crip1 Sept3, Lgals3, Tlr11, Aif1, Dbn1, Btla, Cxcl9, Slc8b1, Ece1, Cst3 Negative: Ccl5, S100a4, Ccnd1, Relb, Ffar2, Siglecg, Dtx1, Ltb, H2-DMb2, Il4i1 Plxnc1, Tbc1d4, Cd7, Cyp4f16, Tmem176b, Traf1, Ppp1r14a, Tmem176a, H2-Eb2, Ryr3 Man1a, Adam23, Lyst, Bcl2a1b, Avpi1, Icosl, Tctex1d2, Arap2, St8sia6, Ube2e2 Computing nearest neighbor graph Computing SNN Computing nearest neighbors Only one graph name supplied, storing nearest-neighbor graph only
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 21646 Number of edges: 732713 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.9031 Number of communities: 22 Elapsed time: 3 seconds
17:33:05 UMAP embedding parameters a = 0.9922 b = 1.112 17:33:06 Commencing smooth kNN distance calibration using 1 thread 17:33:07 Initializing from normalized Laplacian + noise 17:33:08 Commencing optimization for 200 epochs, with 567640 positive edges 17:33:31 Optimization finished 17:33:31 UMAP embedding parameters a = 0.9922 b = 1.112 17:33:31 Read 21646 rows and found 90 numeric columns 17:33:31 Using Annoy for neighbor search, n_neighbors = 20 17:33:31 Building Annoy index with metric = euclidean, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 17:33:33 Writing NN index file to temp file /tmp/RtmpMWeusO/file31487b7b8bd4 17:33:33 Searching Annoy index using 1 thread, search_k = 2000 17:33:37 Annoy recall = 100% 17:33:38 Commencing smooth kNN distance calibration using 1 thread 17:33:39 Initializing from normalized Laplacian + noise 17:33:40 Commencing optimization for 200 epochs, with 682638 positive edges 17:34:05 Optimization finished
reduction <- "rna_umap_nno"
dplot_1 <- DimPlot(so, reduction = reduction, group.by = "rna_snn_res.0.8", label = TRUE) &
theme(aspect.ratio = 1, legend.position = "none")
dplot_2 <- DimPlot(so, reduction = reduction, group.by = "main_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_3 <- DimPlot(so, reduction = reduction, group.by = "fine_labels", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_4 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$tissue, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_5 <- DimPlot(so, reduction = reduction, group.by = "treatment", label = FALSE) &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$treatment, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
fplot_1 <- FeaturePlot(so, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
options(repr.plot.width = 15, repr.plot.height = 10)
dplot <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + plot_layout(ncol = 3)
dplot
sessionInfo()
R version 4.1.0 (2021-05-18) Platform: x86_64-conda-linux-gnu (64-bit) Running under: CentOS Linux 7 (Core) Matrix products: default BLAS/LAPACK: /home/fdeckert/bin/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.15.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] grid stats graphics grDevices utils datasets methods [8] base other attached packages: [1] ComplexHeatmap_2.8.0 viridis_0.6.1 viridisLite_0.4.0 [4] gridExtra_2.3 patchwork_1.1.1 RColorBrewer_1.1-2 [7] ggplot2_3.3.3 dplyr_1.0.6 SeuratObject_4.0.1 [10] Seurat_4.0.2 loaded via a namespace (and not attached): [1] Rtsne_0.15 colorspace_2.0-1 rjson_0.2.20 [4] deldir_0.2-10 ellipsis_0.3.2 ggridges_0.5.3 [7] circlize_0.4.13 IRdisplay_1.0 GlobalOptions_0.1.2 [10] base64enc_0.1-3 clue_0.3-59 spatstat.data_2.1-0 [13] farver_2.1.0 leiden_0.3.8 listenv_0.8.0 [16] ggrepel_0.9.1 RSpectra_0.16-0 fansi_0.5.0 [19] codetools_0.2-18 splines_4.1.0 doParallel_1.0.16 [22] polyclip_1.10-0 IRkernel_1.2 jsonlite_1.7.2 [25] Cairo_1.5-12.2 ica_1.0-2 cluster_2.1.2 [28] png_0.1-7 uwot_0.1.10 shiny_1.6.0 [31] sctransform_0.3.2 spatstat.sparse_2.0-0 compiler_4.1.0 [34] httr_1.4.2 assertthat_0.2.1 Matrix_1.3-4 [37] fastmap_1.1.0 lazyeval_0.2.2 limma_3.48.1 [40] later_1.2.0 htmltools_0.5.1.1 tools_4.1.0 [43] igraph_1.2.6 gtable_0.3.0 glue_1.4.2 [46] RANN_2.6.1 reshape2_1.4.4 Rcpp_1.0.6 [49] scattermore_0.7 vctrs_0.3.8 nlme_3.1-152 [52] iterators_1.0.13 lmtest_0.9-38 stringr_1.4.0 [55] globals_0.14.0 mime_0.10 miniUI_0.1.1.1 [58] lifecycle_1.0.0 irlba_2.3.3 goftest_1.2-2 [61] future_1.21.0 MASS_7.3-54 zoo_1.8-9 [64] scales_1.1.1 spatstat.core_2.1-2 promises_1.2.0.1 [67] spatstat.utils_2.1-0 parallel_4.1.0 reticulate_1.20 [70] pbapply_1.4-3 rpart_4.1-15 stringi_1.6.2 [73] S4Vectors_0.30.0 foreach_1.5.1 BiocGenerics_0.38.0 [76] shape_1.4.6 repr_1.1.3 rlang_0.4.11 [79] pkgconfig_2.0.3 matrixStats_0.59.0 evaluate_0.14 [82] lattice_0.20-44 ROCR_1.0-11 purrr_0.3.4 [85] tensor_1.5 labeling_0.4.2 htmlwidgets_1.5.3 [88] cowplot_1.1.1 tidyselect_1.1.1 parallelly_1.25.0 [91] RcppAnnoy_0.0.18 plyr_1.8.6 magrittr_2.0.1 [94] R6_2.5.0 magick_2.7.2 IRanges_2.26.0 [97] generics_0.1.0 pbdZMQ_0.3-5 DBI_1.1.1 [100] pillar_1.6.1 withr_2.4.2 mgcv_1.8-36 [103] fitdistrplus_1.1-5 survival_3.2-11 abind_1.4-5 [106] tibble_3.1.2 future.apply_1.7.0 crayon_1.4.1 [109] uuid_0.1-4 KernSmooth_2.23-20 utf8_1.2.1 [112] spatstat.geom_2.1-0 plotly_4.9.3 GetoptLong_1.0.5 [115] data.table_1.14.0 digest_0.6.27 xtable_1.8-4 [118] tidyr_1.1.3 httpuv_1.6.1 stats4_4.1.0 [121] munsell_0.5.0